'''
一个天文单位的大小是地球到太阳的平均距离约为1点5亿千米

　　地月距离约为38400千米

　　所以月亮距地球0.000256个天文单位.
太阳、地球和月亮就是一个典型的三体系统，
其中太阳质量为1.989 × 1 0 30 k g 1.989×10^{30}kg1.989×10 
30
 kg，地球质量为5.965 × 1 0 24 k g 5.965×10^{24}kg5.965×10 
24
 kg，月球质量为7.342 ✕ 1 0 22 k g 7.342✕10^{22}kg7.342✕10 
22
 kg，万有引力常数为G = 6.67 × 1 0 − 11 N ⋅ m 2 / k g 2 G=6.67×10^{-11}N·m2/kg^2G=6.67×10 
−11
 N⋅m2/kg 
2
 。地月距离为3.8 × 1 0 8 m 3.8\times10^8m3.8×10 
8
 m；日地距离为1.5 × 1 0 11 m 1.5\times10^{11}m1.5×10 
11
 m；地球公转速度为28.8 k m / s 28.8km/s28.8km/s；月球公转速度为1 k m / s 1km/s1km/s
————————————————
版权声明：本文为CSDN博主「微小冷」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/m0_37816922/article/details/120699335

'''
import time
import numpy as np
from vpython import *

scene.height=880
scene.width=1890



# RE、ME 天文单位
au,G,RE,ME = 1.48e11,6.67e-11,1.48e11,5.965e24
m1=1.989e30
m2=ME
m3=0.000256*ME
r3_init=3.8e8
zoom=500  #一天文单位显示为
zoomTime=3600*24*10  #1秒代表10天
b1 = sphere(pos=vector(0,0,0),radius = (3/40/np.pi*m1)**(1/3),color=color.red,make_trail=True,trail_radius=1)  #定义球的大小和颜色
b2 = sphere(pos=vector(RE,0,0), radius = (0.2/40/np.pi*m1)**(1/3),color=color.blue,make_trail=True,trail_radius=1)  #定义球的大小和颜色
b3 = sphere(pos=vector(r3_init+RE,0,0), radius = (0.2/40/np.pi*m1)**(1/3),color=color.green,make_trail=True,trail_radius=1)  #定义球的大小和颜色
v2_init=28.8e3
v3_init=1e3
b1v=[0,0,0]
b2v=[0,v2_init,0]
b3v=[0,v3_init+v2_init,0]
print(b1v,"    ",b2v,"    ",b3v)
f1=(0,0,0)  # r→g:F1
f2=(0,0,0)  # r→b:F2
f3=(0,0,0)  # g→b:F3
dt0=1/100
dt=dt0*zoomTime
i=0
scene.camera.zoom = 1/RE*zoom/2
scene.camera.pos = vector(0, 0, RE*2)

scene.camera.follow(b1)
while 1:
    # l1 = curve(b1.pos, b2.pos,radius=1)
    # l2 = curve(b1.pos, b3.pos,radius=1)
    # l3 = curve(b2.pos, b3.pos,radius=1)
    time.sleep(dt0)
    
    v1 = (b2.pos.x-b1.pos.x,b2.pos.y-b1.pos.y,b2.pos.z-b1.pos.z)
    v2 = (b3.pos.x - b1.pos.x, b3.pos.y - b1.pos.y, b3.pos.z - b1.pos.z)
    v3 = (b3.pos.x - b2.pos.x, b3.pos.y - b2.pos.y, b3.pos.z - b2.pos.z)
    r1 = (v1[0] ** 2+v1[1]**2+v1[2]**2)**0.5
    r2 = (v2[0] ** 2 + v2[1] ** 2 + v2[2] ** 2) ** 0.5  #原yy
    r3 = (v3[0] ** 2 + v3[1] ** 2 + v3[2] ** 2) ** 0.5

    F1 = G * m1 * m2 / r1/r1
    F2 = G * m1 * m3 / r2/r2
    F3 = G * m2 * m3 / r3/r3

    b1v[0] += (F1 * v1[0] / r1 + F2 * v2[0] / r2) / m1 * dt
    b1v[1] += (F1 * v1[1] / r1 + F2 * v2[1] / r2) / m1 * dt
    b1v[2] += (F1 * v1[2] / r1 + F2 * v2[2] / r2) / m1 * dt
    b2v[0] += (F1 * (-v1[0]) / r1 + F3 * v3[0] / r3) / m2 * dt
    b2v[1] += (F1 * (-v1[1]) / r1 + F3 * v3[1] / r3) / m2 * dt
    b2v[2] += (F1 * (-v1[2]) / r1 + F3 * v3[2] / r3) / m2 * dt
    b3v[0] += (F2 * (-v2[0]) / r2 + F3 * (-v3[0]) / r3) / m3 * dt
    b3v[1] += (F2 * (-v2[1]) / r2 + F3 * (-v3[1]) / r3) / m3 * dt
    b3v[2] += (F2 * (-v2[2]) / r2 + F3 * (-v3[2]) / r3) / m3 * dt
    if i%100==0:
        print(i)
        print(b1v,b2v,b3v)
        print(b1.pos,b2.pos,b3.pos)

    b1.pos.x += b1v[0] * dt
    b1.pos.y += b1v[1] * dt
    b1.pos.z += b1v[2] * dt
    b2.pos.x += b2v[0] * dt
    b2.pos.y += b2v[1] * dt
    b2.pos.z += b2v[2] * dt
    b3.pos.x += b3v[0] * dt
    b3.pos.y += b3v[1] * dt
    b3.pos.z += b3v[2] * dt

    # l1.clear()
    # l2.clear()
    # l3.clear()
    i+=1
    scene.title="b1:"+str((round(b1.pos.x,1),round(b1.pos.y,1),round(b1.pos.z,1)))+"   b2:"+str((round(b2.pos.x,1),round(b2.pos.y,1),round(b2.pos.z,1)))+"   b3:"+str((round(b3.pos.x,1),round(b3.pos.y,1),round(b3.pos.z,1)))
    # if i % 600 == 0:
    #     scene.camera.follow(b1)
    # if i % 600 == 200:
    #     scene.camera.follow(b2)
    # if i % 600 == 400:
    #     scene.camera.follow(b3)
    scene.caption=("Distance R ↔ B ：" +str(round(r1,1)))
    scene.append_to_caption("       Distance R ↔ G ：" + str(round(r2,1)))
    scene.append_to_caption("       Distance B ↔ G ：" + str(round(r3,1)))

